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Abstract 

We have presented some practical consequences on the molecular- 
dynamics simulations arising from the numerical algorithm published re- 
cently in paper pQ. The algorithm is not a finite-difference method and 
therefore it could be complementary to the traditional numerical integrat- 
ing of the motion equations. It consists of two steps. First, an analytic 
form of polynomials in some formal parameter A (we put A = 1 after all) 
is derived, which approximate the solution of the system of differential 
equations under consideration. Next, the numerical values of the derived 
polynomials in the interval, in which the difference between them and 
their truncated part of smaller degree does not exceed a given accuracy 
e, become the numerical solution. The particular examples, which we 
have considered, represent the forced linear and nonlinear oscillator and 
the 2D Lennard- Jones fluid. In the latter case we have restricted to the 
polynomials of the first degree in formal parameter A. 

The computer simulations play very important role in modeling mate- 
rials with unusual properties being contradictictory to our intuition. The 
particular example could be the auxetic materials. In this case, the ac- 
curacy of the applied numerical algorithms as well as various side-effects, 
which might change the physical reality, could become important for the 
properties of the simulated material. 
PACS:31.15.Qg, 02.60.Cb, 02.60.-x 

1 Introduction 

Recently, we have published a numerical algorithm for the Cauchy problem for 
the ordinary differential equations pQ. We showed that it could be much more 

'author: c-mail: mdudek@proton.if.uz.zgora.pl 



1 




Figure 1: The exact solution, x(t) — \ t sin(t) + cos(t), of the forced harmonic 
oscillator, x" + x = cos(t) where £0 = 1, i>o = 0, (the thick line) and three pairs 
of approximating expressions xjsr, (1) and (1') starting at point (Po), (2) and 
(2') starting at point (Pi), (3) and (3') starting at point (P2). We have chosen 
N = 5 in the approximations (1),(2),(3), whereas N = 4 in (1'), (2'), (3')- In 
this example the number of exact digits is equal to 6, e = 10~ 4 . 



accurate, even by few orders of magnitude, than traditional numerical methods 
based on finite differences. In physical applications, the requirement of one 
force evaluation per time step makes that the most often chosen algorithm is 
the Verlet algorithm [US], being the simple third order Taylor predictor method, 
or the equivalent leap-frog algorithm [31 H] • In this case, the possibility to use 
algorithm being much more accurate then Verlet algorithm and as fast as the 
Verlet algorithm makes new perspective for simulating such complex systems as, 
e.g., tetratic phases [5] or auxetics [6]-[8]. Apart from the problem of numerical 
accuracy there is also the possibility of the loss of the time-reversibility in finite- 
difference methods [5], |10) . 

In the following, we discuss our algorithm with respect to integrating the 
motion equations. To this aim we have introduced a few examples of the forced 
linear and nonlinear oscillators and 2D Lennard-Jones fluid. 

2 Short description of the algorithm 

We present the procedure T] of finding an approximate solution of the following 
initial value problem for the second order differential equation of the form: 

x" = f(x,v) + g(t), (1) 

x(t )=x , x'(t ) = v , (2) 
where x' = v, / and g are given functions, and xq, «q are fixed reals. For the 
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Figure 2: The exact solution of the forced harmonic oscillator, x" + x — cos(wi), 
(the thick line), the approximate solution generated by the Velocity- Ver let al- 
gorithm with step size h = 0.01 and the solution generated by our polynomial 
algorithm (of order of A 5 ). In this example e = 10 -2 and oj = 0.1. 

function / we assume that it is sufficiently smooth, so we can write /, using the 
Taylor formula, in some neighborhood of (xq,vq) in the form 

f(x,v) = Z% =0 [(x-x )-^+(v-vo)-^} k f(x ,vo)+o{W(x - *o) 2 + (v- v )Y)- 

(3) 

We introduce a formal real parameter A and instead of the Eqs. UH]) we consider 
the family of problems 

x" = X(f(x,v)+g(t)) (4) 

with the initial data in Eq. @. Next, we seek the approximate solution of 
Eq. (Q in the form 

x N {t) =x + v T + '£% =1 \ k (p k (T), (5) 
where ifk (t) are unknown functions of r = t — to satisfying the condition 

<p k (0) = 4(0) = 0. (6) 

Putting Eq. ([5]) into Eq. ((4]) and next comparing the coefficients of order A fe we get 
the system of differential equations for ip^ , which, together with initial conditions 
Eq.®, determine ifk in the unique way. The differential equations for if/, we 
solve by simple integration. 

To illustrate this procedure we consider the mathematical pendulum problem 
with external force cos(t) 

x" = -x + cos(f), x(to) = xq, x'{t ) = v . (7) 
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For TV = 3 we get 

x 3 = x + v t + E! =1 AVfe(Y), 



(8) 



x'i = EjU^Vfc (r) = -A (so + u r) - S^ =1 A fe +Vfc(r) + Acos(t + t), (9) 

where we substitute to + t for t and the derivatives with respect to t for the 
derivatives with respect to r. Hence, 



p"( T )=- x o-VQT + coa(to + T), <P2(t) = -<pi(r), tp'sir) = -(p 2 (t), (10) 
and after integrating the above equations in the interval [0, r] we obtain 



= cos(i ) - r sin (to) - cos(t + r) - ^— ^ - ^— ^, (11) 



V2(t) = cos(t ) ~ r sin (t ) - cos(t + r) 

r 3 sin(t ) T 5 t; r 4 x 
6 120 24 ' 

93 3 (r) = cos(t ) - t sin(to) - cos(t + t) 
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r 2 cos(to) 
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(12) 



(13) 



We claim that for sufficiently large N and A = 1 the expression xjv(t) is 
a good approximation of the solution of the Eqs. (jl!2p on a small interval of 

f e [t 0l to + <y. 

Practically, for a fixed N we look for the interval of t S [to, to + Si] such that 

\x N (t) -x N -i(t)\ < e, (14) 

where e > is a fixed accuracy. In the above example of mathematical pendulum 
the condition states that \(ps(t — fo)| < e. Next, we repeat our procedure for the 
Eq. (fl]) with the new initial data 

x(t + 6 1 ) = x N (t + 6 1 ), x'(t + 5 1 )=x' N (t + 5 1 ), (15) 

and so on. Fig[T]is a visualization of the updating procedure for the initial data. 
Thus, every time the condition in Eq. (fT4"j) fails at some value of t\ > to, the new 
initial data {to,Xo,i>o} are defined at ti, i.e., to = ti, 2:0 = x(ti), vq — v(t\). 

In many examples it is enough to put N = 3 to get the good approximation 
of the solution. 



4 




algorithm step (Verlet a.) , accuracy (polynomial a.) 



Figure 3: Dependence of the total calculation time of the Velocity- Verlet al- 
gorithm and polynomial method on the chosen value h (in case of the Verlet 
a.) and the given accuracy e (in case of polynomial method). Time units, t v , 
represent calculation time of the Velocit- Verlet algorithm with h = 0.01, 




Figure 4: Two periodic solutions of the forced Duffing oscillator with the same 
parameters a = 0.1 and b = 3.5 and different initial values of xq and vq. The 
attractors have been plotted for t > 400. The parameters of the polynomial 
algorithm: e = 0.6, N = 2. 
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Figure 5: The entire trajectory of the forced Duffing oscillator with the same 
parameters as in Fig. [Hand xq = —1.5 and vq = 0. 

3 Some features of the algorithm 

While performing numerical integrating motion equation one always is fighting 
for the numerical accuracy. In classical finite-difference methods like Verlet 
algorithm, leapfrog algorithm or Runge-Kutta algorithm this is connected with 
the chosen size h of the time step. However, the smaller step size the larger 
cumulated round-off error because more time steps are necessary to cover the 
given time interval. Thus, one should use a numerical method using the smaller 
number of steps (the larger value of h) without loss in numerical accuracy. The 
advantage of our method is already evident in Fig. [2j where three solutions of 
the forced oscillator equation 

x" = -x + cos(wt), x(0) = l, v(0)=0, (16) 

have been plotted, the exact one represented by the equation 

x(t) = - 1 J cos(ujt) - u? cos(i)), (17) 

1 — U! Z 

and two numerical approximations represented by the Velocity- Verlet algorithm 
with the step size h = 0.01 and our polynomial a; at of degree N = 5 in formal 
variable A. In the case of the polynomial method there have been plotted, in the 
figure, only the dots representing the points, where the condition Eq. (|14p fails 
for a given accuracy e = 0.01. They are the only points where the numerical 
round-off errors contribute to the approximate solution. The remaining points 
(in between), which have not been plotted, do not contribute to round-off errors 
cumulation. One always can recalculate them from the exact expression for the 
polynomial representation of x(t). 

The following advantage of our algorithm could be relatively shorter total 
calculation time than in any numerically stable finite-difference method in the 
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Figure 6: Data representing 60ps of the kinetic energy per particle for the 2D 
Lennard- Jones fluid consisting of 500 particles. The chosen step size h = 0.00001 
in the case of the Verlet algorithm whereas e = 0.0001 and N = 1 in the case 
of the polynomial method. 



limit of small values of h. In Fig. [3J we have presented calculation time depen- 
dence of the Velocity- Verlet algorithm on the value of h and our polynomial 
algorithm on the given accuracy e = h. The results in the figure have been 
obtained from the programs calculating deviations of the approximate solutions 
from the exact one. The numerical errors arising from the assumed value of e 
can be by a few orders smaller than in classic finite-difference methods. This 
feature has already been discussed in our paper pQ , where we compared various 
numerical algorithms with respect to their numerical accuracy. 

The next feature of the presented algorithm is that it applies also for strongly 
non-linear motion equations. In particular, in Fig|4j we have presented two 
different attractors of the forced Duffing oscillator (the parameters have been 
taken from the Fig. 2.20 in a book by Holden [IT] ) 

x" + ax' + x 3 = bcos(t), (18) 

with the same values of a = 0.1 and b — 3.5 but different initial conditions. 
In Fig. O there has been presented entire trajectory starting from the initial 
condition and leading to one of the attractors. 

One can use our method also for chaotic solutions of the oscillator. However, 
we do not discuss this possibility in this paper. 

The polynomial X2(t), in the case of the Duffing oscillator, is represented by 
the following formulae: 

x 2 (t) = x + tvq + A [-- t 2 v a - - x r 4 - cos(t + r) b - i s„t 2 

1 3 
~ 2 x l r3 v o + cos(£ ) b - — t 5 vl - t b sin(to)] 
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+A 2 [24 x sin(t + r)bv Q + 72 sin(i + r) r b Ug + ^ x 3 , t 6 vf 
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— 12 xo cos(io) r 6 ug — 6 sin(to) a + ^ x o 7 b sin(to ) 

— 108 cos(^o) 6 Vg — 3xg r6 sin(to) — 12 x r cos(t + r) bvo 
-^xl cos(to) t 2 b + ~ Xg 5 t 3 a 

+ 108 cos(i + t) 6u 2 + 36 r 6 up sin(to) - 3xg cos(i + r) 6 



+- t 5 -r 4 
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+ 40 X ° ^ W ° ~~ 2 COS ^ ^ r4 6w o + 2 r2 6 sin (*o) « 
53 

-24 x 6u sin(to) + y^jjj Xg r 7 «g + sin(i + r) 6 a 

9 1 

+ 10 1-5 ^ U " sin ( io ) + - r 3 w a 2 + t 4 6 v Q sm(t ) 

2 3 

- x r 5 Up a - 18 t 2 cos(t + r) b v% + — r 9 v% 



1 3 

+ - Xg T 4 Vg a + 3 Xg cos(tg) b — cos(tg) t 6 a + - s t 5 JJq] 



(19) 



and in this case the accepted approximated solution should satisfy the inequality 
^(t)! < e for a given e, where v?2(V) is the coefficient of A 2 (we set the formal 
parameter A = 1 after all). 

In our previous paper [1] we have shown that our method could be used also 
for molecular-dynamics simulations of large number of particles. To this end, we 
have simulated barometric formula in the case of the ideal gas of 1000 molecules 
in the gravitational field and the gas was contacted Nose-Hoover thermostat 

H2, [TJ. 

In all mentioned by us cases, till now, the series expansion of the force 
(Eq. ((3|)) consisted of a finite number of terms. The question arises, could 
the method be extended to a more general case, where the number of terms is 
infinite? In order to show the possibility we have considered 2D Lennard-Jones 
fluid represented by a system of N particles interacting with Lennard-Jones 
potential energy, 



U( nj ) = 4e 



(-)--(-) 



(20) 



Then, the force experienced by the particle i from another particle j being a 
distance away is repesented by the following formula 



dU(r, 
dru 



24e 
a 



2(-) i3 -e) 7 



(21) 



In this case the series expansion in the neighborhood of [xq, Vq) (see Eq.[3]) leads 
to an infinite number of terms including the powers of \J (xoi — xoj ) 2 + (yoi — Uoj) 2 - 
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In the case of the approximating polynomials of the order of A 1 the numerical 
algorithm is equivalent to the Velocity- Verlet algorithm and it is represented by 
the following set of equations: 

► 1 ^ > T T 2 

n{t)=r^ + v^r + —{¥ l +F lV -) — , (22) 

^)=vm + -t(fI + F^?-), (23) 
m 2 

where 

Fi = (r$-r]$)\Ft\, (24) 

FiV = (vlo - Vjo)\Fi\ ■ (25) 

and fio and ■u^o are the initial location and velocity of the particle i, respectively. 
The accuracy control parameter e should satisfy the condition 

l^ + ^)yl<- (26) 

The generalization of the algorithm to the case of the polynomial approximation 
of the order A 2 becomes much more complex and it is not presented in this 
paper. However, already the results obtained in the linear approximation (in 
formal parameter A) become promising. In Fig. [6l there has been plotted the 
kinetic energy (per particle) of 500 particles representing 2D Lennard-Jones 
fluid versus time in the case of the Velocity- Verlet algorithm and our polynomial 
approximation (Eqs. I22ll23p . linear in A. In this case, the total time tp used by 
the polynomial method was of the same order as the total time ty of the Verlet 
method (tp — 1.5 ty). In order to preserve the given numerical accuracy e the 
polynomial algorithm was runnining according to the following steps: 

(i) start with t = ho, where ho = 0.001 

(ii) if Eq. (|26|) fails then change the value of r by some factor, e.g. r = r/10 

(iii) calculate the values of the polynomials rt and vt 

(iv) goto (i). 

In the considered example, the time intervals r used during the entire simulation 
run were distributed as follows: 



T = 


10" 3 


- 0% 


T = 


io- 4 


= 49% 


r = 


10~ 5 


= 48.7% 


T = 


10" 6 


= 2%. 



(27) 



9 



The total calculation time strongly depends on the value of e and the higher 
order of the approximating polynomial (in the formal variable A) makes possible 
larger values of r. 

4 Conclusions 

We have discussed possible numerical advantages of integrating motion equa- 
tions with the help of the recently published algorithm for solving the initial- 
value problem for the ordinary differential equations [T] . Contrary to the tradi- 
tional finite-difference methods, representing truncated series expantion of the 
solution of the equation of motion under consideration, the method is not dis- 
crete in time. This makes possible that for the large class of problems in physics 
the algorithm could be faster and more accurate than traditional finite-difference 
schemes. The particular example of the 2D Lennard-Jones fluid, which has been 
discussed above, suggests that the method could be applied to the many-body 
problems 
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